clear;clc;
A=[0.98 0.1 0;0.02 0.7 0.05;0 0.2 0.95];
T=1000000;
states=zeros(1,T);
[b,v]=eig(A);

%(a)
%State C1 is represented by 1,state C2 by 2 and state O by 3,respectively
states(1)=2;
for t=1:T-1
   states(t+1)=myfun(states(t),A);
end
figure;
hist(states,6);
title('the distribution of states');
xlim([0 4]);
set(gca,'XTick',[1 2 3]);

%(b)
a1=length(find(states==1));
a2=length(find(states==2));
a3=length(find(states==3));
A1=1;A2=a2/a1;A3=a3/a1;
B1=1;B2=b(2,2)/b(1,2);B3=b(3,2)/b(1,2);
figure;plot([1,2,3],[A1 A2 A3],'o',[1 2 3],[B1 B2 B3],'+');
axis([0 4 0 1.5]);
legend('eigenvector','states','Location','best');
title('comparasion between states and dominant eigenvector with relative value');

%(c)
%State "open" is represented by -0.5,state "close" by 0.5,respectively
rstates=zeros(1,T);
for i=1:T
    if states(i)<3
        rstates(i)=-0.5;
    else
        rstates(i)=0.5;
    end;
end;
figure;
hist(rstates,5);
title('the distribution of rstates');
xlim([-1 1]);
set(gca,'XTick',[-0.5 0.5]);

%(d)
t = 800; 
open=zeros(1,t);
closed=zeros(1,t);
 m=1;
 max=0;
for k=2:T
    if (rstates(k)==rstates(k-1))
        m=m+1;
    else  
       if (rstates(k-1) == 0.5)
           closed(m)=closed(m)+1;
       else
           open(m)=open(m)+1;
       end;
       if (m>max) max=m;
       end
       m=1;
    end
end;

figure;bar(open);
figure;bar(closed);

%(e)
